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ABSTRACT 

. An efficient simulation algorithm for random sequential adsorption of discs is implemented. 

. By dividing the surface into small squares, depositions are attempted only on squares that are 

actually possible. A crucial part is a method to identify the available squares. A precise value of 
i-C . the jamming coverage is obtained: d{oo) w 0.547069. The asymptotic law 9{t) w 0{oo) — ci^^/^ is 

["t I . verified to a high accuracy. The pair correlation function at the jamming state is analyzed. 

l> 

1. Introduction 

Random sequential adsorption (RSA) is an irreversible process which models depositions of large 
■ molecules (e.g. proteins) on surfaces.^ In RSA large extended objects are deposited on a flat surface one by 

\^ ' one, with locations chosen at random. The deposited objects are fixed on the surface and exclude certain 

o ■ 

' hard discs on a two-dimensional continuum.-^ ^ Depositions of squares, elongated objects like ellipses 

and needles, mixtures of different shapes, ^^'^^ as well as lattice models^'^^'^°~^'^ were also studied by 
many authors. Experiments on RSA lagged behind theoretical studies, but beginning to catch up.^^'^^ 

RSA exhibits several interesting features not present in equilibrium systems. First of all, there exists a 
jamming state, which has the highest possible particle density smaller than the close-packed density of the 
pH ■ corresponding equilibrium system. At the jamming state, the distances between objects are such that no 

O ■ other object can be added to the surface. The jamming density is approached via a power law on continuum 

and an exponential law on discrete lattices. At the jamming state, the pair correlation function diverges 
logarithmically when approaching the distance of contact. 
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^ • The coverage (or density) at large time t takes the form 



0(t) « 0(oo) - (1) 

where c > is some constant, the exponent df depends on the geometry of the objects. The most interesting 
quantities are the value d{(X)) and the scaling exponent df. Analytic arguments^^'^^ and numerical work^'^ 
support the result that df = d, for hyperspheres in d dimensions. In two dimensions the jamming limit is 
approaches by c/v^ for discs. This result is referred as Feder's law in the literature. For anisotropic objects 
like ellipses and unoriented squares, ^"'^^'^^ the jamming coverage is approached even slower, d/ = 3. It is 
argued that d/ is the number of degrees of freedom of the object in general.^" On the other hand, for the 
deposition of squares with fixed orientation (oriented squares), the convergence law is^^'i^'^v i^ji^^^ 

In this paper, we present an efficient method of the RSA simulation and report the results of extensive 
runs on a cluster of workstations. 
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2. Standard and Event-Driven Simulation of RSA 

The standard RSA model of hard discs^"^ is defined as follows: We begin with an empty, flat surface of 
L X L at time t = 0. A point (x, y) is drawn from a uniform probability distribution over the entire surface 
area. A deposition of a disc of diameter a = 1 with center at {x, y) is attempted. If this disc does not 
overlap with previously deposited discs, i.e., if the distances between the current point {x,y) and the centers 
of all other discs on the surface are greater than a, the deposition is accepted, otherwise it is rejected. Each 
attempt increases the time by 1/LP'. The deposition trial is repeated until no more deposition is possible, at 
which point, the jamming state is reached. The coverage 6{t) is the fraction of the area covered by the discs 
at time t. 

The standard computer simulation of RSA follows closely the description of the model. The checking 
for overlap can be done quickly with the help of grids. ^ The total L x L surface area is partitioned into 
regions of squares of size 1 x 1. A disc at {x,y) belongs to the unit square at ([a;J, [j/J), where the notation 
[.xj represents the floor or integer part of x. An array of linked lists of disc coordinates is used to represent 
the deposited discs. With this data structure, it is suSice just to check the discs located at the center and 
eight surround squares for overlaps, see Fig. 1. 

The above simple method has a severe shortcoming that the study of late-stage process is very time- 
consuming. When the jamming configuration is about to be reached, most of the area is blocked. Only 
tiny disconnected regions are available for further deposition. Since we choose the sites at random with 
a uniform probability distribution, the available sites are harder to find. To overcome this difficulty, we 
incorporated two important ideas: (1) divide the surface into small squares and make deposition attempts 
only on squares that are not completely blocked; (2) systematically reduce the sizes of the small squares 
after some number of attempts and re-evaluate the availability of the squares. The first strategy is similar 
to that for the simulation of Ising model at low temperatures,^* as well as methods of RSA.^^'^^ The second 
strategy appears new and very effective. 

In this more elaborate method, we make attempts only on squares which are potentially possible for 
depositions. Thus the core part is an algorithm which identifies correctly whether a square is available for 
deposition. Since the diameter of the discs is a, each disc excludes an area of circle of radius a for further 
deposition (Fig. 1). If some area is completely covered by a union of the exclusion zones of discs, that area 
is not available for deposition. We begin by classiiying all the squares of a x a with a = 1 as available or 
unavailable squares. The available squares are put on a list. Deposition trials are taken only on the squares 
in the list, i.e., a square in the list is chosen at random, deposition is attempted with a uniform probability 
over the square. After certain number of trials (we used 5 x 10^), we rebuild the list, now with squares of 
size (a/2) x (a/2), checking only those squares on the old list. Each old square is subdivided into 4 smaller 
squares. This shrinking of the basic squares helps to locate even extremely small available regions. They 
will be hit with very small probability if the area is always 1x1. Because of integer overflow, the size of the 
squares is not allowed to go arbitrarily small, but the shrinking is stopped at a = 2~^^. Thus the smallest 
square has an area of2~^° w 10~^. 

For the standard RSA algorithm, each trial increases the time t by 1/A, where A = is the area of the 
total surface. In the even-driven algorithm, where only the potentially successful depositions are tried, the 
clock must tick faster in order to compensate for not making attempt in the completely unsuccessful area. 
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The time increment for each trial on the available squares is a random variable/^ 

In^ 



-4 



ln(l- AM). 



+ 1 , (2) 



where is the total area on which we try our deposition, it is equal to times the number of available 
squares; and ^ is a uniformly distributed random number between and 1. This interpretation of time 
interval makes the two algorithms equivalent. 

The above result can be understood in the following way: the total area is divided into Ag and A — Ag. 

If we try to deposit the discs uniformly on the whole area A, the probability that a trial landed on Ag is 
p — As/ A. If the first attempt fails to hit Ag, a second attempt landed on Ag has the probability p(l — p). 
In general, the probability that the first « — 1 trials landed in the unavailable area A — As and the i-th trial 
landed on Ag is 

P,=pil-py-\ z=l,2, 3,---. (3) 

In the event-driven simulation, from first trial in the area As to a second trial in the same area, the time 
has elapsed by At = i/A, since i — I unsuccessful attempts had been made in the area A — Ag, where i is 
distributed according to Eq. 3. This exponential distribution is easily generated by the logarithmic function, 
given Eq. 2. 

3. Identiflcation of Fully Covered Squares 

An important piece of code of our RSA simulation program is the identification of the squares which 
are fully covered by the excluded area of discs. For these squares we are sure that depositions on them will 
be unsuccessful, and thus they'll not be on the list of candidates for trials. 

The algorithm that we have devised is as follows. Written as a C programming language function, it 
returns a nonzero value if the square is fully covered and a value otherwise. Part (1) to (5) is executed in 
the order given. 

(1) Find relevant discs to the current square. A disc is relevant if the square overlaps with the excluded 
region of the disc. Note that a disc has a diameter a and its excluded region is a concentric circle of 
radius a. 

(2) If all of the four vertices of the square are inside a single-disc excluded region, then it is already fully 
covered (Fig. 2a). Function returns with value 1. 

(3) For a full coverage, each vertex of the square must be at least in one of the excluded region of the 
relevant discs. If any one of the vertices is not covered at all by any disc, then the square is not fully 
covered (Fig. 2b). Function returns with value 0. 

(4) At this point, at least two discs are involved if the function is not returned. Mapping the edges of the 
square to a one-dimensional line between and 4, we find out all the line segments which are covered 
by the excluded region of discs. If the four edges of the square arc completely covered by the discs, and 
we have exactly two discs, then the square is fully covered (Fig. 2c), function returns with value 2. If 
there are segments not covered by discs, then we know that the square is not fully covered (Fig. 2d). 
Function returns with value 0. 

(5) If the function is not returned, we know that there are at least three discs. Even if all the edges are 
covered, the square can still contain uncovered area if three or more discs are involved. Three or more 
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discs can form holes. To make the last check, the coordinates of the intersection points of circles (with 
radius a) are calculated. In order for the interior of the square to be covered, all the intersection points 
of any pair of discs which lie inside or on the square edges must be covered by a third disc. Return with 
the number of relevant discs if the above condition is satisfied (Fig. 2c); return value if not (Fig. 2f). 

The order of various checks is arranged in such a way so that the most frequent sitiiations arc checked first. 

Even though the worse case computational complicity goes as ri^, where n is the number of relevant discs, 

a definite conclusion can be made typically much earlier. 

4. Simulation Results 

We used mainly a system of linear size L — 1000. This appears to be the largest system considered in 
RSA simulation. With this size various data occupied 24 megabytes. Periodic boundary condition is used. 
This is crucial to eliminate an obvious edge effect of order 1/L. It is known that finite-size eSect in RSA is 
extremely small^^ due to rapid damping of two-point correlations. The actual form of the finite-size effect is 
not known. In any case, a system of 1000 x 1000 is more than adequate and finite-size effect can be neglected. 

For time t <8we used the standard method. This method is simpler, thus faster in the initial deposition 
process. The use of the event-driven algorithm at this stage would require much more memory. The event- 
driven method is used for t > 8. In the beginning, the size of small squares (for classification purpose) is 
1x1, which is the same as the underline grids for book-keeping of adjacency of discs. The squares are 
subdivided into four squares after every 0.5L^ trials. This helps to keep the number of available squares less 
than 0.3L^, even though the area of the squares is smaller. 

A simulation terminates if there is no more available square, or if t > 2^^ (a complete run). On a total 
area of 1000 x 1000, the jamming state is reached typically at a time t « 2*^ with large fiuctuation. In general, 
the jamming state is obtained on the time scale of t « i^. In each complete process of deposition of discs 
on a finite lattice, the coverage varies discretely. The last disc deposited on the surface causes a variation 
in the coverage by an amount of order Comparing this variation with the result of an infinitely large 

system or with the result of average over many runs, which is of order t~^^^, we conclude that the major 
deposition process is finished when these two numbers are of the same order of magnitude, given t ^ . 

Each complete run took 22 minutes on an SGI Indigo workstation. The standard simulation method 
ran at the speed of about one step {At = 1, attempts) per minute, independent of the time t. If we were 
determined to use this method for all t, a complete run would have taken 2^^ minutes (« 8 million years) to 
finish. 

Floating-point calculations are carried out in double precision. We used two random number gener- 
ators in the program, a 64-bit linear congruential generator^^ and a 48-bit one implemented on many C 
programming language environment on workstations. Both of them take the form 

Xn+i = axn + cmod m, (4) 

with a = 6 364136 223 846 793005, c = 0, m = 2^^ for the 64-bit random number and a = 25 214903 917, 
c = 11, m = 2*^ for the 48-bit random number. 

The coverage data are taken at t = 2'^, with k = 0, 1, 2, 49. The pair correlation function is 
calculated from the last configuration. The final results are averages over 43 427 runs. 
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To check the convergence law of the coverage, we plot in Fig. 3 log2[^(2t) — 6{t)\ vs. log2f. 

Assuming a power-law dependence of the form in Eq. 1, we have 




where cq = log2[c(l — 2~^/'^})\. By a linear least-squares fit, we find that df = 2.0008 ± 0.0016 and 
c = 0.300 ± 0.002, confirming the theoretical expectation [df = 2) to high accuracy. 

To get the jamming coverage, we can simply take the value at t = 2^^. The actual difference between 
t = 2^^ and t — > oo is in fact much smaller than the statistical errors. The estimate for the jamming coverage 
is thus 



This result agrees with previous work,^"'*'^'^^ 9{oo) « 0.547, but is two to three orders of magnitude more 
accurate. 

The pair correlation function of the final configuration (at t = 2^^) is presented in Fig. 4. The pair 
correlation function g(r) is defined as the average number of discs per unit area located at a distance r away, 
given that there is a disc at the origin, normalized by the jamming density so that (7(00) = 1. A striking 
feature is the divergence at r — > 1. Generalizing an exact result in one dimension, Pomeau^^ and Swendsen^^ 
showed that 



The logarithmic divergence is nicely confirmed from our numerical data (Fig. 4, insert). A least-squares fit 
with data ln(r - 1) < -6 gives a = -1.535 ± 0.002 and b = -2.13 ± 0.02. 

5. Conclusion 

We proposed a very efficient simulation algorithm for the random sequential adsorption of discs. With 
this new method and extensive runs, accurate jamming coverage is obtained. The Feder's law for the large 
time asymptotic coverage and logarithmic divergence of the pair correlation are confirmed to a high accuracy. 
The method can be easily generalized to RSA of other geometric shapes with fixed orientation. It appears 
complicated to apply to problems where the objects have a rotational degree of freedom, e.g., deposition of 
ellipses. 



9{oo)= 0.547 069 ± 0.000 000 7. 



(6) 



g{r) « aln(r — 1) + 6. 



(7) 
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Figure Captions 

Fig. 1. RSA configuration with five discs deposited on surface. Each disc has an exclusion zone (black and 
grey region). The grids are used to locate the neighboring discs. 

Fig. 2. Various geometric relations between the square and excluded regions of discs. In situation (a), (c), 
or (e) the square is fully covered by one, two, or more than two discs; (b), (d), or (f) fails to cover the 
square. 

Fig. 3. The difference in coverage, 9{2t) — 0{t), vs. time t, plotted in double logarithmic (to the base 2) 
scale. The data points are simply connected by straight lines. 

Fig. 4. Pair correlation function g{r) at the jamming state. Insert is a portion near r — > 1. 
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